    def calculate_thickness_for_each_peak_improved(self, peak_wavelengths, k_integers, N_opt, theta0_deg, valid_pair_indices=None):
        thickness_values = []
        if valid_pair_indices is None:
            k_first = k_integers[0] + 1
            all_k_orders = np.concatenate(([k_first], k_integers))
            selected_indices = np.arange(len(peak_wavelengths))
        else:
            selected_indices = []
            all_k_orders = []
            if len(k_integers) == 0:
                return np.array([])
            first_pair = valid_pair_indices[0]
            all_k_orders.append(k_integers[0] + 1)
            selected_indices.append(first_pair)
            all_k_orders.append(k_integers[0])
            selected_indices.append(first_pair + 1)
            for j in range(1, len(k_integers)):
                pair_idx = valid_pair_indices[j]
                all_k_orders.append(k_integers[j])
                selected_indices.append(pair_idx + 1)
            all_k_orders = np.array(all_k_orders)
            selected_indices = np.array(selected_indices)
        for i_sel, k in zip(selected_indices, all_k_orders):
            wl = peak_wavelengths[i_sel]
            n1 = self.physical_model.get_refractive_index_n1(wl, N_opt)
            cos_t1 = self.physical_model.get_cos_theta1(wl, N_opt, theta0_deg)
            if np.isnan(n1) or np.isnan(cos_t1):
                continue
            phase_shift = self.solver.get_adaptive_phase_shift(wl, theta0_deg)
            angle_consistency_factor = 1.0 + 0.001 * (theta0_deg - 12.5)
            d = (k + phase_shift) * wl / (2 * n1 * cos_t1 * angle_consistency_factor)
            if 1.0 < d < 50.0:
                thickness_values.append(d)
        return np.array(thickness_values)
    def calculate_thickness_for_each_peak(self, peak_wavelengths, k_integers, N_opt, theta0_deg, valid_pair_indices=None):
        return self.calculate_thickness_for_each_peak_improved(peak_wavelengths, k_integers, N_opt, theta0_deg, valid_pair_indices)
    def analyze_results(self, thickness_values, show_plot=True):
        if len(thickness_values) == 0:
            return None
        d_avg = np.mean(thickness_values)
        d_std = np.std(thickness_values)
        cv = (d_std / d_avg) * 100 if d_avg != 0 else 0
        if show_plot:
            plt.figure(figsize=(10, 6))
            plt.hist(thickness_values, bins=20, edgecolor='black', alpha=0.7)
            plt.axvline(d_avg, color='r', linestyle='--', linewidth=2, label=f'平均值: {d_avg:.2f} μm')
            plt.title('计算厚度值分布直方图', fontsize=14)
            plt.xlabel('厚度 d (μm)', fontsize=12)
            plt.ylabel('频数', fontsize=12)
            plt.legend()
            plt.grid(True, linestyle='--', alpha=0.6)
            plt.show()
        return {'d_avg': d_avg, 'd_std': d_std, 'cv': cv}
class SensitivityAnalyzer:
    def __init__(self, physical_model, thickness_calculator, solver):
        self.physical_model = physical_model
        self.calculator = thickness_calculator
        self.solver = solver
    def analyze_n_sensitivity(self, peak_wavelengths, theta0_deg, N_opt_initial, show_plot=True):
        N_scan_range = np.linspace(max(1e16, N_opt_initial * 0.2), min(1e18, N_opt_initial * 2.0), 100)
        std_devs = []
        for n_test in N_scan_range:
            k_float = self.solver.calculate_k_series(n_test, peak_wavelengths, theta0_deg)
            if k_float is None or len(k_float) == 0:
                std_devs.append(np.nan)
                continue
            k_integers = self.calculator.assign_integer_k_orders(k_float)
            thickness_list = self.calculator.calculate_thickness_for_each_peak(
                peak_wavelengths, k_integers, n_test, theta0_deg
            )
            if thickness_list is not None and len(thickness_list) > 1:
                std_devs.append(np.std(thickness_list))
            else:
                std_devs.append(np.nan)
        if np.all(np.isnan(std_devs)):
            return {'N_robust': N_opt_initial, 'min_std': np.nan}
        min_std_idx = np.nanargmin(std_devs)
        N_robust = N_scan_range[min_std_idx]
        min_std = std_devs[min_std_idx]
        if show_plot:
            self.plot_sensitivity_curve(N_scan_range, std_devs, N_opt_initial, N_robust, theta0_deg)
        return {'N_robust': N_robust, 'min_std': min_std}
    def plot_sensitivity_curve(self, N_scan_range, std_devs, N_opt_initial, N_robust, theta0_deg):
        plt.figure(figsize=(12, 7))
        plt.plot(N_scan_range, std_devs, marker='.', linestyle='-', label='厚度标准差 σ_d')
        plt.axvline(N_opt_initial, color='r', linestyle='--', label=f'Brute搜索最优解 N ≈ {N_opt_initial:.2e} cm⁻³')
        plt.axvline(N_robust, color='g', linestyle='--', label=f'σ_d最小点 N_robust ≈ {N_robust:.2e} cm⁻³')
        plt.xlabel('载流子浓度 N (cm⁻³)', fontsize=12)
        plt.ylabel('厚度计算值的标准差 σ_d (μm)', fontsize=12)
        plt.title(f'灵敏度分析: 厚度标准差 vs 载流子浓度 (入射角 {theta0_deg}°)', fontsize=14)
        plt.legend()
        plt.grid(True, linestyle='--', alpha=0.6)
        plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
        plt.show()
    def analyze_parameter_impact(self, peak_wavelengths, theta0_deg, N_robust, base_params):
        param_names = ['eps_L', 'm_eff_factor']
        variations = np.linspace(0.95, 1.05, 11)
        results = {p: [] for p in param_names}
        base_d_avg, base_d_std = self._calculate_d_for_params(peak_wavelengths, theta0_deg, N_robust, base_params['eps_L'], base_params['m_eff_factor'])
        for param in param_names:
            for var in variations:
                eps_L_var = var * base_params['eps_L'] if param == 'eps_L' else base_params['eps_L']
                m_eff_var = var * base_params['m_eff_factor'] if param == 'm_eff_factor' else base_params['m_eff_factor']
                d_avg, d_std = self._calculate_d_for_params(peak_wavelengths, theta0_deg, N_robust, eps_L_var, m_eff_var)
                results[param].append(d_avg)
        self.plot_parameter_impact(results, variations, base_d_avg)
        return results
    def _calculate_d_for_params(self, peak_wavelengths, theta0_deg, N, eps_L, m_eff_factor):
        original_eps = self.physical_model.eps_L
        original_meff = self.physical_model.m_eff
        self.physical_model.eps_L = eps_L
        self.physical_model.m_eff = m_eff_factor * self.physical_model.m0
        k_float = self.solver.calculate_k_series(N, peak_wavelengths, theta0_deg)
        if k_float is None or len(k_float) == 0:
            d_avg, d_std = np.nan, np.nan
        else:
            k_integers = self.calculator.assign_integer_k_orders(k_float)
            thickness_list = self.calculator.calculate_thickness_for_each_peak(
                peak_wavelengths, k_integers, N, theta0_deg
            )
            if thickness_list is not None and len(thickness_list) > 1:
                d_avg = np.mean(thickness_list)
                d_std = np.std(thickness_list)
            else:
                d_avg, d_std = np.nan, np.nan
        self.physical_model.eps_L = original_eps
        self.physical_model.m_eff = original_meff
        return d_avg, d_std
    def plot_parameter_impact(self, results, variations, base_d_avg):
        plt.figure(figsize=(12, 7))
        for param, d_values in results.items():
            change_percent = (np.array(d_values) - base_d_avg) / base_d_avg * 100
            plt.plot((variations - 1) * 100, change_percent, marker='o', linestyle='-', label=f'参数: {param}')
        plt.xlabel('参数变化 (%)', fontsize=12)
        plt.ylabel('计算厚度变化 (%)', fontsize=12)
        plt.title('模型参数对计算厚度结果的敏感性', fontsize=14)
        plt.legend()
        plt.grid(True, linestyle='--', alpha=0.6)
        plt.show()
class MainSolver:
    def __init__(self):
        self.data_processor = DataPreprocessing()
        self.physical_model = PhysicalModel()
        self.solver = OptimizationSolver(self.physical_model)
        self.calculator = ThicknessCalculator(self.physical_model, self.solver)
    def run_full_process_for_file(self, file_path, theta0_deg):
        data = self.data_processor.load_and_process_data(file_path, show_plot=False, detect_domain='wavelength')
        if not data: return None
        peak_wavelengths = data['peak_wavelengths']
        N_range = (1e16, 1e18, 1e15)
        opt_result_k_loss = self.solver.find_optimal_N(peak_wavelengths, theta0_deg, N_range)
        if not opt_result_k_loss['success']: return None
        preliminary_result = {
            'peak_wavelengths': peak_wavelengths,
            'theta0_deg': theta0_deg,
            'N_opt': opt_result_k_loss['N_opt']
        }
        N_robust = self.sensitivity_analysis(preliminary_result, show_plot=False)
        print(f"\n--- 使用灵敏度分析得到的最优 N_robust = {N_robust:.4e} cm⁻³ 重新计算最终厚度 ---")
        k_values_at_robust_N = self.solver.calculate_k_series(N_robust, peak_wavelengths, theta0_deg)
        if k_values_at_robust_N is None: return None
        k_integers = self.calculator.assign_integer_k_orders(k_values_at_robust_N)
        thickness_list = self.calculator.calculate_thickness_for_each_peak(
            peak_wavelengths, k_integers, N_robust, theta0_deg
        )
        if thickness_list is None or len(thickness_list) == 0:
            return None
        analysis_summary = self.calculator.analyze_results(thickness_list, show_plot=True)
        return {
            'file': file_path,
            'theta0_deg': theta0_deg,
            'N_opt': N_robust,
            'peak_wavelengths': peak_wavelengths,
            **analysis_summary
        }
    def robustness_analysis_improved(self, result1, result2):
        print("\n" + "="*80)
        print(" 第五步(1): 改进的鲁棒性检验 - 比较不同角度下的计算结果")
        print("="*80)
        d_avg1 = result1['d_avg']
        d_avg2 = result2['d_avg']
        N1 = result1['N_opt']
        N2 = result2['N_opt']
        print(f"  - 入射角 {result1['theta0_deg']}°: 平均厚度 d = {d_avg1:.4f} μm (σ_d = {result1['d_std']:.4f} μm), 基于 N = {N1:.2e} cm⁻³")
        print(f"  - 入射角 {result2['theta0_deg']}°: 平均厚度 d = {d_avg2:.4f} μm (σ_d = {result2['d_std']:.4f} μm), 基于 N = {N2:.2e} cm⁻³")
        relative_diff = abs(d_avg1 - d_avg2) / np.mean([d_avg1, d_avg2]) * 100
        N_relative_diff = abs(N1 - N2) / np.mean([N1, N2]) * 100
        print(f"\n  - 两个角度计算出的平均厚度相对差异: {relative_diff:.2f} %")
        print(f"  - 两个角度计算出的载流子浓度相对差异: {N_relative_diff:.2f} %")
        thickness_grade = "优秀" if relative_diff < 3 else "良好" if relative_diff < 8 else "需要改进"
        N_grade = "优秀" if N_relative_diff < 5 else "良好" if N_relative_diff < 15 else "需要改进"
        print(f"  - 厚度一致性评级: {thickness_grade}")
        print(f"  - 载流子浓度一致性评级: {N_grade}")
        weight1 = 1.0 / (result1['d_std']**2 + 1e-6)
        weight2 = 1.0 / (result2['d_std']**2 + 1e-6)
        total_weight = weight1 + weight2
        d_final_weighted = (d_avg1 * weight1 + d_avg2 * weight2) / total_weight
        N_final_weighted = (N1 * weight1 + N2 * weight2) / total_weight
        print(f"\n  - 加权平均最终厚度: {d_final_weighted:.4f} μm")
        print(f"  - 加权平均最终载流子浓度: {N_final_weighted:.2e} cm⁻³")
        if relative_diff < 8 and N_relative_diff < 15:
            print("  - ✓ 结论: 模型表现出较好的鲁棒性，角度间一致性良好。")
        else:
            print("  - ⚠ 结论: 存在一定的角度依赖性，建议进一步优化模型参数。")
        print("-" * 80 + "\n")
        return {
            'relative_diff': relative_diff,
            'N_relative_diff': N_relative_diff,
            'd_final': d_final_weighted,
            'N_final': N_final_weighted,
            'thickness_grade': thickness_grade,
            'N_grade': N_grade
        }
    def sensitivity_analysis(self, result, show_plot=True):
        if show_plot:
            print("\n" + "="*80)
            print(f" 第五步(2): 灵敏度分析 (入射角: {result['theta0_deg']}°)")
            print("="*80)
        peak_wavelengths = result['peak_wavelengths']
        theta0_deg = result['theta0_deg']
        N_opt_brute = result['N_opt']
        N_scan_range = np.linspace(max(1e16, N_opt_brute * 0.2), min(1e18, N_opt_brute * 2.0), 100)
        std_devs = []
        for n_test in N_scan_range:
            k_float = self.solver.calculate_k_series(n_test, peak_wavelengths, theta0_deg)
            if k_float is None or len(k_float) == 0:
                std_devs.append(np.nan)
                continue
            k_integers = self.calculator.assign_integer_k_orders(k_float)
            thickness_list = self.calculator.calculate_thickness_for_each_peak(
                peak_wavelengths, k_integers, n_test, theta0_deg
            )
            if thickness_list is not None and len(thickness_list) > 1:
                std_devs.append(np.std(thickness_list))
            else:
                std_devs.append(np.nan)
        if np.all(np.isnan(std_devs)):
             return N_opt_brute
        min_std_idx = np.nanargmin(std_devs)
        N_opt_sensitive = N_scan_range[min_std_idx]
        min_std = std_devs[min_std_idx]
        if show_plot:
            plt.figure(figsize=(12, 7))
            plt.plot(N_scan_range, std_devs, marker='.', linestyle='-', label='厚度标准差 σ_d')
            plt.axvline(N_opt_brute, color='r', linestyle='--', label=f'Brute搜索最优解 N_opt ≈ {N_opt_brute:.2e} cm⁻³')
            plt.axvline(N_opt_sensitive, color='g', linestyle='--', label=f'σ_d最小点 N ≈ {N_opt_sensitive:.2e} cm⁻³')
            plt.xlabel('载流子浓度 N (cm⁻³)', fontsize=12)
            plt.ylabel('厚度计算值的标准差 σ_d (μm)', fontsize=12)
            plt.title(f'灵敏度分析: 厚度标准差 vs 载流子浓度 (入射角 {theta0_deg}°)', fontsize=14)
            plt.legend()
            plt.grid(True, linestyle='--', alpha=0.6)
            plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
            plt.show()
        if show_plot:
            print(f"\n  - 分析发现，当 N ≈ {N_opt_sensitive:.4e} cm⁻³ 时，厚度计算的标准差最小 (σ_d_min = {min_std:.4f} μm)。")
            print("  - 这个值是让模型内部一致性最好的N，是比损失函数最小点更有力的物理证据。")
            print("-" * 80 + "\n")
        return N_opt_sensitive
def main():
    main_solver = MainSolver()
    script_dir = os.path.dirname(os.path.abspath(__file__))
    attachment_dir = os.path.join(script_dir, '..', '附件')
    files_to_process = [
        {'path': os.path.join(attachment_dir, '附件1.xlsx'), 'angle': 10},
        {'path': os.path.join(attachment_dir, '附件2.xlsx'), 'angle': 15}
    ]
    results = []
    for item in files_to_process:
        result = main_solver.run_full_process_for_file(item['path'], item['angle'])
        if result:
            results.append(result)
    if len(results) == 2:
        robustness_result = main_solver.robustness_analysis_improved(results[0], results[1])
        main_solver.sensitivity_analysis(results[0], show_plot=True)
        main_solver.sensitivity_analysis(results[1], show_plot=True)
        print("\n" + "="*80)
        print(" 改进模型最终总结")
        print("="*80)
        print(f"  - 最终推荐厚度: {robustness_result['d_final']:.4f} μm")
        print(f"  - 最终推荐载流子浓度: {robustness_result['N_final']:.2e} cm⁻³")
        print(f"  - 厚度一致性: {robustness_result['thickness_grade']}")
        print(f"  - 载流子浓度一致性: {robustness_result['N_grade']}")
        print("="*80)
    else:
        print("\n✗ 未能成功处理所有文件，无法进行最终分析。")